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pvj Abstract As the title suggests, the purpose of this chapter is to review the current 

, status of numerical simulations of black hole accretion disks. This chapter focuses 

O . exclusively on global simulations of the accretion process within a few tens of grav- 

itational radii of the black hole. Most of the simulations discussed are performed 
using general relativistic magnetohydrodynamic (MHD) schemes, although some 
mention is made of Newtonian radiation MHD simulations and smoothed particle 
hydrodynamics. The goal is to convey some of the exciting work that has been 
going on in the past few years and provide some speculation on future directions. 
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1 Introduction 



X 
6 

"*^ Going hand-in-glove with analytic models of accretion disks, discussed in Chap- 
el ter 2.1, are direct numerical simulations. Although analytic theories have been 
'~^ extremely successful at explaining many general observational properties of black 
^^ hole accretion disks, numerical simulations have become an indispensable tool in 
►^ advancing this field. They allow one to explore the full, non-linear evolution of 
y—{ accretion disks from a first-principles perspective. Because numerical simulations 
■^^ can be tuned to a variety of parameters, they serve as a sort of "laboratory" for 
l/~J astrophysics. 

'^ The last decade has been an exciting time for black hole accretion disk simula- 

^^ tions, as the fidelity has become sufficient to make genuine comparisons between 

^^ them and real observations. The prospects are also good that within the next 

C^ decade, we will be able to include the full physics (gravity -|- hydrodynamics + 
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Fig. 1 Model light c urve at 0.4 mm (solid) and accretio n rate (dotted) at the inner boundary 
of a simulation from [McKinney and Blandford| | |2009[ | . Bot h quantities are sca led to their 
maximum value. These are compared to data from Sgr A* in [Dexter et al.|||2010[l. 



magnetic fields + radiation) within these simulations, which will yield complete 
and accurate numerical representations of the accretion process. In the rest of this 
chapter I will review some of the most recent highlights from this field. 



2 Matching simulations with observations 



One of the most exciting recent trends has been a concerted effort by various col- 
laborations to make direct connections between very sophisticated simulations and 
observations. Of course, observers have been clamoring for this sort of comparison 
for years! 

Perhaps the first serious attempt at this was presented in [Schnittman et aI7 
(2006'). Schnittman produced a simulation similar to those in De Villiers et al. 
(2003) and coupled it with a ray-tracing and radiative transfer code to produce 
"images" of what the simulated disk would look like to a distant observer. By 
creating images from many time dumps in the simulation, Schnittman was able to 
create light curves, which were then analyzed for variability properties much the 
way real light curves are. 

Following that same prescription, a number of groups have now presented 
results coupling general relativistic MHD (GRMHD) simulations with radiative 
modeling and ray-tracing codes (e.g. Noble et al.||2007 Moscibrodzka et al.||2009 



Dexter et al. 2009[ |2010[). More recent models have even included polarization 
measurements ( Shcherbakov et al.|2012 1 . This approach is most applicable to very 
low-luminosity systems, such as Sgr A* and M87. A sample light curve for Sgr A* 
covering a 16-hour window is shown in Figure [l| In the case of M87, modeling has 
focused on accounting for the prominent jet in that system (Moscibrodzka et al. 
20 11 [ [Dexter et al.||2012[ ). 
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Fig. 2 Synt hetic broadband spectrum created from one of the simulations presented in 

et al.| ||2012[|. The pin k points represent a compilation of Sgr A* observations. Figure from 
Drappeau et al!]||2013[l. 



Along with modeling light curves and variability, this approach can also be 
used to create synthetic broadband spectra from simulations (e.g. 'Moscibrodzka| 



et al. 2009 



Drappeau et al.^_2013J , which can be compared with modern multi- 
wavelength observing campaigns (see Chapter 3.1). This is very useful for con- 
necting different components of the spectra to different regions of the simulation 
domain. For example, Figure [2] shows that the sub-mm bump in Sgr A* is well 
represented by emission from relatively cool, high-density gas orbiting close to the 
black hole, while the X-ray emission seems to come from Comptonization by very 
hot electrons in the highly magnetized regions of the black hole magnetosphere or 
base of the jet. 



3 Thermodynamics of simulations 



As important as the radiative modeling of simulations described in Section [2] has 
been, its application is very limited. This is because, in most cases, the radia- 
tive modeling has been done after the fact; it was not included in the simulations 
themselves. Therefore, the gas in the accretion disk was not allowed to respond 
thermodynamically to the cooling. This calls into question how much the structure 
obtained from the simulation reflects the true structure of the disk. Fortunately, 
various groups are beginning to work on treating the thermodynamics of accre- 
tion disks within the numerical simulations with greater fidelity. Thus far, two 
approaches have principally been explored: 1) ad hoc cooling prescriptions used to 
artificially create optically thick, geometrically thin disks and 2) fully self-consistent 
treatments of radiative cooling for optically thin, geometrically thick disks. We review 
each of these in the next 2 sections. 
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Fig. 3 Various fluxes as functions of radius for a numerical Novikov-Thorne disk simulation. 
Top: Mass accretion rate. Second panel: Accreted specific angular momentum. Solid line is 
simulation data; dashed line gives Novikov-Thorne solution; dotted line is ISCO value. Note 
that the specific angular momentum does not drop significantly inside the ISCO. Third panel: 
The "nominal" efficiency, which is the total loss of specific energy from the fluid. Bottom panel: 
Speciflc magnetic flux. The near constancy of this quantity inside t he ISCO is an indic ation 
that magnetic stresses are not significant in this region. Figure from|Penna et ar]l|2010}. 



3.1 Geometrically thin disks 

For the ad hoc cooling prescription, cooling is assumed to equal heating (approx- 
imately) everywhere locally in the disk. Since this is the same assumption as is 
made in the Shakura-Sunyaev ( Shakura and Sunyaev|[l973 1 and Novikov-Thorne 
(Novikov and Thorne 19731 disk models, this approach has proven quite useful 



in testing the key assumptions inherent in these models (e.g Shafee et al 



2008 



Noble et al.||2009[ |Penna et al.||2010[ [Noble et al.||2010[ ). In particular, these simu- 



lations have been useful for testing the assumption that the stress within the disk 
goes to zero at the innermost stable circular orbit (ISCO). A corollary to this is 
that the specific angular momentum of the gas must remain constant at its ISCO 
value inside this radius. Both of these effects have been confirmed in simulations 
of sufficiently thin disks ( Penna et al.||2010 ) , as shown in Figure [s] 



3.2 Self-consistent radiative cooling of optically thin disks 



Another approach to treating the thermodynamics of accretion disks has been 
to include physical radiative cooling processes directly within the simulations. So 
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Fig. 4 Comparison of sample spectra generated for Sgr A* from numerical simulations at 
three different accretion rates: 10~® (black), 10~* (blue), and W~''^ Mq yr~^ (red). For each 
accretion rate, two simulations are shown, one that includes cooling self-consistently (model 
names ending in "C") and one that does not. T he spectra begin to diverge noticeably at 
M Ri 10-^Mq yr-^ Figure from|Dibi et al.||2012l. 



far there has been very hmited work done on this for opticaUy thick disks, but 



an optically-thin treatment was introduced in Fragile and Meier (2009). Similar 
to the after-the-fact radiative modeling described in Section l2j the optically-thin 
requirement restricts the applicability of this approach to relatively low luminosity 
systems, such as the Quiescent and Low/Hard states of black hole X-ray binaries. 
Recently this approach has been applied to Sgr A* (Drappeau et al. 20131, 



which turns out to be right on the boundary between where after-the-fact ra- 
diative modeling breaks down and a self-consistent treatment becomes necessary 
( Dibi et al.||2012 ). Figure [4] illustrates that this transition occurs right around an 
accretion rate of M « 1O~^M0 yr~^ for Sgr A*. 



4 Magnetic field topology 

Another area where a lot of interesting new results have come out is in the study 
of how magnetic field topology and strength affect black hole accretion. 



4.1 Jet power 

Although there is now convincing evidence that the Blandford-Znajek mechanism 



Komis- 



(Blandford and Znajek 1977) works as predicted in powering jets (e.g 
sarov||2001 McKinney and Gamniie||2004 ), one lingering question is still how the 
accretion process supplies the required poloidal fiux onto the black hole. Simu- 
lations have demonstrated that such field can, in many cases, be generated self- 



consistently within MRI-unstable disks (e.g. De Villiers et al 



2005 



Hawley and 
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Fig. 5 Time-ave raged siiell integ r als of data from unbound outflows as a function of radius 
for 3 models from [Beckwith et al.| | [2008| : a dipole magnetic field model (black), a quadrupole 
model (cyan), and a toroidal field model (magenta) . Shown are mass outflow rate M (top left), 
magnetic field strength ||fe'^|| (top right), electromagnetic energy flux |rj'"|EM (bottom left), and 
angular momentum flux ITTIbm (bottom right). Dashed lines show ±1 standard deviation from 
the average. 



Krolik 2006 McKinney and Gammie 2004 McKinney 20061. However, this is 



strongly dependent on the initial magnetic field topology, as shown in |Beckwith| 
(20081. Figure [5] nicely illustrates that when there is no net poloidal mag- 



et al. 



netic flux threading the inner disk, the magnetically-driven jet can be 2 orders of 
magnitude less energetic than when there is. At this time it is unclear what the 
"natural" field topology would be, or even if there is one. 



4.2 Magnetically arrested accretion 



Although strong poloidal magnetic fields are useful for driving powerful jets, they 
can also create interesting feedback affects on an accretion disk. In the case where 
a black hole is able to accumulate field with a consistent net fiux for an extended 
period of time, it is possible for the amassed field to eventually "arrest" the ac- 
cretion process ( Narayan et al.|2003 ). An example of an arrested state is shown in 
Figure |6) 

In a two-dimensional simulation where plasma with a constant net flux is fed 
in from the outer boundary, a limit-cycle behavior can set in, where the mass 
accretion rate varies by many orders of magnitude between the arrested and non- 
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Fig. 6 Distributions of density in the meridional plane at different simulation times, showing 
a magnetically-arrested state {leff) and a non-arrested state (right). Botto m: Snapshot of 
magnetic field lines at the same simulation times. Figure from |Igumenshchev| pOOS^ . 



arrested states. Figure [7] provides an example of the resulting mass accretion his- 
tory. It is straightforward to show that the interval, At, between each non-arrested 
phase in this scenario grows with time according to 



At ~ '' r 



(1) 



where Vr is the radial infall velocity of the gas, Bz is the strength of the magnetic 
field, and p is the density of the gas. 

In three-dimensions, the magnetic fields are no longer able to perfectly arrest 
the in- falling gas because of a "magnetic" Ray leigh- Taylor effect. Basically, as low 
density, highly magnetized gas tries to support higher density gas in a gravitational 
potential, it becomes unstable to an interchange of the low- and high-density 
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Fig. 7 Evolution of mass accretion rate and magnetic fluxes in 2D (axisymmetric) simulation 
of magnetically-arrested accretion. Accretion into the black hole begins at t Si 1.3. Starting 
from t ■^ 1.4, a pattern o f cyclic accretion develops (seen as a sequence of spikes). Figure 
from |Igumenshchev| ||2008[l . 



materials. Indeed, such a magnetic Ray leigh- Taylor effect has been seen in recent 



(2012) 



simulations by Igumenshchev (20081; Tchekhovskoy et al. (2011); McKinney et al. 



The results of Tchekhovskoy et al. 



(20111 are important for another reason. 
These were the first simulations to demonstrate a jet efhciency 77 = (M — E)/{M) 
greater than unity. Since the efhciency measures the amount of energy extracted 
by the jet, normalized by the amount of energy made available via accretion, a 
value ri > 1 indicates more energy is being extracted than is being supplied by 
accretion. This is only possible if some other source of energy is being tapped - in 
this case the rotational energy of the black hole. This was the first demonstration 
that a Blandford-Znajek ( Blandford and Znajek|1977 l process must be at work in 
driving these simulated jets. 



5 Tilted disks 



There is observational evidence that several black-hole X-ray binaries (BHBs), 
e.g. GRO J1655-40 (|Orosz and Bail^pOT] ), V4641 Sgr ( jMiher et al.][2OT2| and 



GX 339-4 (Miller et al. 20091, and active galactic nuclei (AGN), e.g. NGC3079 



(|Kondratkoetal.|2005[ ), NGC 1068 ( |Caproni et al.|2006[ ), and NGC 4258 ( [Caproni 



et al.|2007 l, may have accretion disks that are tilted with respect to the symmetry 
plane of their central black hole spacetimes. There are also compelling theoretical 
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Fig. 8 (Panels a-d): the top and bottom rows show, respectively, equatorial (z = 0) and 
meridional {y = 0) snapshots of the gas density of a magnetically arrested flow in 3D. Black 
lines show field lines in the image plane. (Panel e): time evolution of the mass accretion rate. 
(Panel f): time evolution of the large-scale magnetic flux thre ading the BH horizon . (Pane l g): 
time evolution of the energy outflow efficiency rj. Figure from |Tchekhovskoy et aL](|201l|. 



arguments that many black hole accretion disks should be tilted (Fragile et al. 2001 
Maccarone|2002 )■ This applies to both stellar mass black holes, which can become 
tilted through asymmetric supernovae kicks ( Fragos et al.|2010 ) or binary captures 
and will remain tilted throughout their accretion histories, and to supermassive 
black holes in galactic centers, which will likely be tilted for some period of time 
after every major merger event ( Kinney et al.|p000 |. 

Close to the black hole, tilted disks may align with the symmetry plane of the 



black hole, either through the Bardeen-Petterson effect (Bardeen and Petterson 



1975 1 in geometrically thin disks or through the magneto-spin alignment effect 



(McKinney et al. 2013) in the case of geometrically thick, magnetically-choked 
accretion (McKinney et al. 2012). However, for weakly magnetized, moderately 



thick disks {H/r > 0.1), no alignment is observed (Fragile et al. 2007). In such 
cases, there are many observational consequences to consider. 



5.1 Tilted disks and spin 



Chapter 4.3 of this book discusses the two primary methods for estimating the 
spins of black holes: continuum-fitting and reflection-line modeling. Both rely on 
an assumed monotonic relation between the inner edge of the accretion disk (as- 
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Fig. 9 Plot of the effective inner radius rin of simulated untilted {circles) and tilted (diamonds) 
accretion disks as a function of black-hole spin using a surface dens ity measure E{ri^) = 
i^max/3e. The solid Une is the ISCO radius. Figure from [Fragile! | |2009| . 



sumed to coincide with the radius of the ISCO) and black hole spin, a. This is 
because what both methods actually measure is the effective inner radius of the 
accretion disk, rin. One problem with this is that it has been shown ( Fragile|2009l 

I I I , , , ■ 1 , 1 T 1 1 , r 11 1 , '• 1 1 



Dexter and Fragile"2011) that tilted disks do not follow such a monotonic behav- 



ior, at least not for disks that are not exceptionally geometrically thin. Figure [9] 
shows an example of the difference between how rjn depends on a for untilted and 
tilted simulated disks. Similar behavior has been confirmed using both dynamical 
([Fragile 2009) and radiative (Dexter and Fragile 2011) measures of rin. The im- 
plication is that spin can only be reliably inferred in cases where the inclination 
of the inner accretion disk can be independently determined, such as by modeling 
jet kinematics (Steiner and McClintock]|2012). 



5.2 Tilted disks and Sgr A* spectral fitting 



For geometrically thin, Shakura-Sunyaev type accretion disks, the Bardeen-Petterson 



effect (Bardeen and Petterson 1975) may allow the inner region of the accretion 



disk to align with the symmetry plane of the black hole, perhaps alleviating con- 
cerns about measuring a, at least for systems in the proper state ("Soft" or "Ther- 
mally dominant") and luminosity range L < 0.3Z/Edd, where I/Edd is the Eddington 
luminosity. Extremely low luminosity systems, though, such as Sgr A*, do not ex- 
perience Bardeen-Petterson alignment. Further, for a system like Sgr A* that is 
presumed to be fed by winds from massive stars orbiting in the galactic center, 
there is no reason to expect the accretion flow to be aligned with the black hole 
spin axis. Therefore, a tilted conflguration should be expected. In light of this. 



Dexter and Fragile (20121 presented an initial comparison of the effect of tilt on 



spectral fitting of Sgr A*. Figure 10 gives one illustration of how important this ef- 
fect is; it shows how the probability density distribution of four observables change 
if one simply accounts for the two extra degrees of freedom introduced by even a 
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Fig. 10 Normalized probability distributions as a function of black hole spin (top left), sky 
orientation (top right), inclination (bottom left), and accr etion rate (bo t tom r ight) for untilted 
(red) and tilte d (black) simulations. Figure adapted from |Dexter et ar] l |2010[ | and [Dexter and] 
|Fragile| ( |20T2l l. 



modestly tilted disk. The take away point should be clear - ignoring tilt artificially 
constrains these fit parameters! 



5.3 Tilted disks and strong shocks 



One remarkable outcome of considering tilt in fitting the spectral data for Sgr 
A* is that tilted simulations seem able to naturally resolve a problem that had 
plagued earlier studies. Spectra produced from untilted simulations of Sgr A* have 
always yielded a deficit of fiux in the near-infrared compared to what is observed. 
For untilted simulations, this can only be rectified by invoking additional spec- 
tral components beyond those that naturally arise from the simulations. Tilted 
simulations, though, produce a sufficient population of hot electrons, without any 
additional assumptions, to produce the observed near-infrared flux (see comparison 



in Figure 11 ). They are able to do this because of another unique feature of tilted 



disks: the presence of standing shocks near the line-of- nodes at small radii ( Fragile 
and Blaesl|2008 1 . These shocks are a result of epicyclic driving due to unbalanced 
pressure gradients in tilted disks leading to a crowding of orbits near their respec- 
tive apocenters. Figure [T2| shows the orientation of these shocks in relation to the 
rest of the inner accretion flow. 
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Fig. 11 Spectra from untilted (left) and tilted (right) simulations. Symbols represent Sgr A* 
data. In both cases the spectra are fit to the green sub-min points. In the tilted simulations, 
multiple electron populations, some heated by shocks associate with the tilt, are present and 
can naturally produce the observed NIR emission, which i s underproduced by ^ 2 o rders of 
magnitude in comparable untilted simulations. Figure from [Dexter and Fragile] (|2012^. 




isodensity surface 



Fig. 12 Three-dimensional contours of density (semi-transparent gray) and Lagrangian spe- 
cific entropy generat ion rate in arbitrary units (solid gray), indicating the shocks. Figure from 
[Henisey et al.| | |2012| |. 



5.4 Tilted disks - GRMHD vs. SPH 



A worthwhile future direction to pursue in this area would be a robust comparison 
of tilted disk simulations using both GRMHD and smoothed-particle hydrodynam- 
ics (SPH) numerical methods. The GRMHD simulations (e.g. Fragile et al. 2007} 



2009 



McKinney et al.|2013 1 enjoy the advantage of being "first principles" calcu- 



lations, since they include all of the relevant physics, whereas the SPH simulations 
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Fig. 13 Gas density overlaid with velocity vectors from three different radiation MHD sim- 
ulations of a black hole accretion disk, probing luminosities of L ~ 10 ~^ (left), 10^''(center), 
and ILEdd (right). Figure adapted from|Ohsuga and Mineshige|||2011[|. 



(e.g. Nelson and Papaloizou 2000| Lodato and Pringle 2007 Lodato and Price 



2010 



[Nixon and King 2012) enjoy the advantage of being more computationally 
efficient, though they make certain assumptions about the form of the "viscosity" 
in the disk. Thus far, the GRMHD and SPH communities have proceeded sepa- 
rately in their studies of tilted accretion disks, and it has yet to be demonstrated 
that the two methods yield equivalent results. This would seem to be a relatively 
straightforward and important thing to check. 



6 Future direction - radiation MHD 



A few years ago, it might have been very ambitious to claim that researchers 
would soon be able to perform global radiation MHD simulations of black hole 
accretion disks, yet a lot has happened over that time, so that now it is no longer 
a prediction but a reality. In the realm of Newtonian simulations, a marvelous 



study was published by Ohsuga and Mineshige (2011), showing global (though 
two-dimensional) radiation MHD simulations of accretion onto a black hole in 
three different accretion regimes: L <C -f/Eddj L < L^^^, and L > L^dd- The 
remarkably different behavior of the disk in each of the simulations (illustrated in 



Figure 13 ) is testament to how rich this field promises to be as more groups join 



this line of research. The specifics of this work are discussed more in Chapter 5.3. 
The other big thing to happen (mostly) within the past year is that a number 
of groups have now tackled, for the first time, the challenge of developing codes for 
relativistic radiation MHD in black hole environments (Farris et al. 2008; Zanotti 



et al.l2011[|Roedig et al.|2012[|FVagile et al.|20T2| |Sg.dowski et al. 2013; Takahashi 
et al.| 20131. So far none of these groups have gotten to the point of simulating 
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Fig. 14 Logarithm of the optical depth (left) and ra diative flux (right) of a black hole accreting 
from a wind passing from left to right. Figure from [Zanotti et al.|||201l|l. 



accretion disks in the way Ohsuga and Mineshige (2011 1 did (they are still mostly 
at the stage of code tests and simple one- and two-dimensional problems), but 
with so many groups joining the chase, one can surely expect rapid progress in the 
near future. One result of some astrophysical interest is the study of Bondi-Hoyle 



(wind) accretion onto a black hole, including the effects of radiation (Zanotti et al. 
2011[ [Roedig et al.||2012i ) (see Figure [l4|. 
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